A Personalized Statistical Approach for Detecting Changes in Gait Cycles Due to Fatigue

Department of Mechanical & Aerospace Engineering, University at Buffalo

Department of Mechanical & Aerospace Engineering, University at Buffalo

Department of Industrial and Systems Engineering, University at Buffalo

Farmer School of Business, Miami University

Farmer School of Business, Miami University

---

1 R Setup and Required Packages

In this project, the open-source R programming language is used to model/monitor the changes in gait cycle patterns due to fatigue in simulated manual material handling task experiments. R is maintained by an international team of developers who make the language available at The Comprehensive R Archive Network. Readers interested in reusing our code and reproducing our results should have R installed locally on their machines. R can be installed on a number of different operating systems (see Windows, Mac, and Linux for the installation instructions for these systems). We also recommend using the RStudio interface for R. The reader can download RStudio for free by following the instructions at the link. For non-R users, we recommend the Hands-on Programming with R for a brief overview of the software’s functionality. Hereafter, we assume that the reader has an introductory understanding of the R programming language.

In the code chunk below, we load the packages used to support our analysis. Note that the code of this and any of the code chunks can be hidden by clicking on the ‘Hide’ button to facilitate the navigation. The reader can hide all code and/or download the Rmd file associated with this document by clicking on the Code button on the top right corner of this document. Our input and output files can also be accessed/ downloaded from saebragani/Gait_EWMA.

# check if packages are not installed; if yes, install missing packages
packages = c("tidyverse", "magrittr", # typical data analysis packages
             "MALDIquant", # match closest points between two vectors
             "foreach", "doParallel", "parallel", # packages for parallelization
             "R.matlab", # to read mat files
             "fda.usc", "MFHD", "qcc") # for depth and EWMA calculations
newPackages = packages[!(packages %in% installed.packages()[,"Package"])]
if(length(newPackages) > 0) install.packages(newPackages)

# using the library command to load all packages; invisible used to avoid printing all packages and dependencies used
invisible(lapply(packages, library, character.only = TRUE))

source("./Functions.R") # custom built and modified

---

2 Gait Acceleration Data and Rational Subgroups

Fourteen participants were equipped with an IMU on their right ankle and performed a three hour manual material handling (MMH) task. The participants loaded weighted boxes on a dolly, pushed the dolly on a set path of 80 m, and unloaded the boxes at the starting location to complete one walking cycle. The IMU signals were collected during the 3-hour sessions. We converted the raw IMU acceleration signals from the local to the global reference frame and removed the gravity. The acceleration signals were then transformed back to the local reference frame. Sagittal acceleration, lateral acceleration, and acceleration magnitude signals were then calculated for further analysis. We used the vertical acceleration component in the global reference frame to segment the gait cycles in order to isolate individual gait cycles. Segmented gait acceleration profiles along with the experimental time stamps were stored in mat files in Matlab.

In the following code chunk we load the segmented acceleration mat files. The \(1^{st}\) 10 minutes of the data were considered as warm up period and thus excluded form the analysis. The acceleration profiles of the gait cycles during each walking cycle were grouped into a rational subgroup. We stored the start and end times of the walking cycles in csv files and loaded and used them in the following code chunk. The acceleration profiles within each start and end times to subgroup the gait cycles.

for (id in setdiff(1:15, 13)) {
  #################### Read start and end time of the subgroups
  video <- read.csv(paste0(file="../Data/csvFiles/Sub", id, ".csv"))

  video_leave <- video$Leaves
  video_leave <- video_leave[!is.na(video_leave)]
  video_leave <- video_leave - 600
  video_leave <- video_leave[video_leave > 0]

  video_enter <- video$Enters
  video_enter <- video_enter[!is.na(video_enter)]
  video_enter <- video_enter - 600
  video_enter <- video_enter[video_enter > 0]

  ################################################### Read from mat files
  raw1 <- readMat(paste0("../Data/matFiles/Subject", id, "_aZ_seg.mat"))

  raw2 <- raw1$gait

  num_rows <- length(raw2)/5

  aM <- list("vector")
  aS <- list("vector")
  aL <- list("vector")
  exp_time <- c()
  for (i in 1:num_rows) {
    aM[[i]] <- raw2[[1*num_rows + i]][[1]][1,]
    aS[[i]] <- raw2[[2*num_rows + i]][[1]][1,]
    aL[[i]] <- raw2[[3*num_rows + i]][[1]][1,]
    exp_time[i] <- raw2[[i+4*num_rows]][[1]][1,1]
  }

  ###################################################### Create subgroups
  if (video_leave[1] > video_enter[1]){
    video_leave <- c(exp_time[1], video_leave)
  }

  if (tail(video_leave, 1) > tail(video_enter, 1)){
    video_enter <- c(video_enter, tail(exp_time, 1))
  }

  conf <- 0.05 * mean(video_enter - video_leave) # 5% confidence for each subgroup

  # Match video subroup times against the experimental times from mat files
  leave_index <- match.closest(video_leave + conf, exp_time)
  enter_index <- match.closest(video_enter - conf, exp_time)

  aM_list <- list()
  aS_list <- list()
  aL_list <- list()
  t_list <- list()
  for (i in 1:length(leave_index)) {
    aM_list[[i]] <- aM[leave_index[i]:enter_index[i]]
    aS_list[[i]] <- aS[leave_index[i]:enter_index[i]]
    aL_list[[i]] <- aL[leave_index[i]:enter_index[i]]
    t_list[[i]] <- exp_time[leave_index[i]:enter_index[i]]
  }

  len <- c()
  for (i in 1:length(aM_list)) {
    len <- c(len, lengths(aM_list[[i]]))
  }

  cut_len <- quantile(len, 0.02)

  for (i in 1:length(aM_list)) {
    remove_ind <- which(lengths(aM_list[[i]]) < cut_len)
    remove_ind <- c(remove_ind, 3000)
    aM_list[[i]] <- aM_list[[i]][-remove_ind]
    aS_list[[i]] <- aS_list[[i]][-remove_ind]
    aL_list[[i]] <- aL_list[[i]][-remove_ind]
    t_list[[i]] <- t_list[[i]][-remove_ind]
  }

  batchSize <- lengths(aM_list)

  ###############################################
  aM_list_all <- do.call(c, aM_list)
  aS_list_all <- do.call(c, aS_list)
  aL_list_all <- do.call(c, aL_list)
  t_all <- do.call(c, t_list)
  ###############################################
  aM_mat <- do.call(rbind, aM_list_all)[,1:cut_len]
  aS_mat <- do.call(rbind, aS_list_all)[,1:cut_len]
  aL_mat <- do.call(rbind, aL_list_all)[,1:cut_len]

  ##################
  assign(paste0("sub", id, "_aMSLT"),
         list(aMag=aM_mat, aSag=aS_mat, aLat=aL_mat, time=t_all, batchSize=batchSize))
}

########################################### Save
save(sub1_aMSLT, sub2_aMSLT, sub3_aMSLT, sub4_aMSLT, sub5_aMSLT, sub6_aMSLT,
     sub7_aMSLT, sub8_aMSLT, sub9_aMSLT, sub10_aMSLT, sub11_aMSLT, sub12_aMSLT,
     sub14_aMSLT, sub15_aMSLT,

     file="../Data/rData/subGs_for_depth.Rdata")

3 Depth Calculation

The depth measure investigates how deep a data point is with respect to a sample. In univariate data the median represents the deepest point in a cloud of points. In functional data analysis (FDA) depth measure is commonly used for classification and clustering of curves. In multivariate data where a number of discretized curves in a certain number of points (usually time points) are investigated, the depth measures the centrality or outlyingness of each curve. However, for multivariate data there is no obvious simple method for depth measurement. Multiple depth measures have been proposed in literature Cuevas et al. (2007). These depth values are defined based on median Fraiman & Muniz (2001), Tukey depth, random projections of the curves and their projections.

In the current work we have exploited the depth measure to investigate the centrality of the gait cycle acceleration profiles in order to find their deviation from the baseline due to fatigue. After investigating multiple depth measures Lamooki et al. (n.d.), we have narrowed them down to mode depth @Cuevas et al. (2007) and Multivariate functional halfspace depth (MFHD) Hubert et al. (2012).

3.1 Mode Depth

In the following code chunk we have calculated the mode depth for each gait cycle acceleration profile using the fda.usc package in R Febrero Bande & Oviedo de la Fuente (2012). The mode depth is calculated for the acceleration magnitude by using the \(1^{st}\) 500 gait cycles as the baseline. The depth of the remaining gait cycle acceleration profiles are calculated against this baseline.

load(file="../Data/rData/subGs_for_depth.Rdata")

for (id in setdiff(1:15, 13)) {
  
  aMag <- get(paste0("sub", id, "_aMSLT"))$aMag
  tExp <- get(paste0("sub", id, "_aMSLT"))$time
  
  ################################## Mode Depth
  inControl <- fdata(aMag[1:500,])
  online <- fdata(aMag[501:nrow(aMag),])
  
  inControlDepth <- {depth.mode(inControl)}$dep
  onlineDepth <- {depth.mode(fdataobj=online, fdataori=inControl)}$dep
  
  modeDepth <- c(inControlDepth, onlineDepth)
  
  assign(paste0("sub", id, "_mode_mag"),
         list(modeDepth=modeDepth, tExp=tExp))
}

save(sub1_mode_mag, sub2_mode_mag, sub3_mode_mag, sub4_mode_mag, sub5_mode_mag, sub6_mode_mag,
     sub7_mode_mag, sub8_mode_mag, sub9_mode_mag, sub10_mode_mag, sub11_mode_mag, sub12_mode_mag, sub14_mode_mag, sub15_mode_mag,
     file="../Data/rData/mode_mag.Rdata")

3.2 Visualize the Mode Depth

The Mode depth for every gait cycle is visualized here where the red data points represent the baseline profiles.

load(file="../Data/rData/mode_mag.Rdata")

for (id in setdiff(1:15, 13)) {
  depthMode <- get(paste0("sub", id, "_mode_mag"))$modeDepth
  tMode <- get(paste0("sub", id, "_mode_mag"))$tExp

  cat("###", paste0("Subject", id), "{-}",'\n')
  plot(tMode, depthMode, pch=16, cex=1,col=c(rep("red", 500), rep("black", length(depthMode)-500)))
  legend("topright", legend=c("Baseline data", "New data"), col=c("red", "black"), pch=16)
  cat('\n \n')

}

Subject1

Subject2

Subject3

Subject4

Subject5

Subject6

Subject7

Subject8

Subject9

Subject10

Subject11

Subject12

Subject14

Subject15

3.3 Multivariate Functional Halfspace Depth (MFHD)

The MFHD depth is introduced to provide ordering of curves from the center outwards for multivariate curves Hubert et al. (2012). A multivariate curve contains a multivariate vector at each discrete time point. The method is particularly used on multivariate curves that are constructed by adding univariate curves and their derivatives. However, in the current work we have implemented the method on multivariate curves constructed by sagittal and lateral gait acceleration profiles. The \(1^{st}\) 500 gait profiles were considered to be the baseline data where the participant was walking fatigue-free. Next, each of the remaining gait acceleration profile was added to the baseline profiles and their MFHD depth value were calculated.

The following code chunk includes the code that we have used to calculate th MFHD depth using the MFHD function in R. Since the MFHD depth computations are significantly time consuming, we have performed the computations in parallel on a supercomputing cluster with 40 cores and set the eval option to FALSE here. The output of this code chunk is added to the folder Data/rData.

load(file="../Data/rData/subGs_for_depth.Rdata")

cores = detectCores() - 2 # to give the server some breathing Room
cl = makePSOCKcluster(cores)
registerDoParallel(cl)

for (id in setdiff(1:15, 13)) {
  
  aSag <- get(paste0("sub", id, "_aMSLT"))$aSag
  aLat <- get(paste0("sub", id, "_aMSLT"))$aLat
  tExp <- get(paste0("sub", id, "_aMSLT"))$time
  
  ################################## Mode Depth
  inControlSag <- aSag[1:500,]
  inControlLat <- aLat[1:500,]
  inControlDepth <- {MFHD(y1=inControlSag, y2=inControlLat, alpha=0.125, Beta=0.5)}$MFHDdepth[1,]
  
  ###################### foreach
  onlineDepth <- c()
  end <- nrow(aSag)
  onlineDepth <- foreach(i=501:end, .packages = c('fda.usc', 'tidyverse', 'MFHD'),
                         .combine='c') %dopar% {
                           temp1 <- aSag %>% .[i,]
                           temp2 <- aLat %>% .[i,]
                           
                           append1 <- rbind(inControlSag, temp1)
                           append2 <- rbind(inControlLat, temp2)
                           
                           temp3 <- MFHD(y1=append1, y2=append2, alpha=0.125, Beta=0.5)
                           
                           temp4 <- temp3$MFHDdepth %>% .[1,501]
                           temp4
                         }
  
  MFHDdepth <- c(inControlDepth, onlineDepth)
  
  assign(paste0("sub", id, "_MFHD"),
         list(MFHDdepth=MFHDdepth, tExp=tExp))
}

save(sub1_MFHD, sub2_MFHD, sub3_MFHD, sub4_MFHD, sub5_MFHD, sub6_MFHD, sub7_MFHD, sub8_MFHD,
     sub9_MFHD, sub10_MFHD, sub11_MFHD, sub12_MFHD, sub14_MFHD, sub15_MFHD,
     file="../Data/rData/MFHD.Rdata")

3.4 Visualize the MFHD Depth

Below we have visualized the MFHD depth for the gait cycles where the red data points represent the baseline gaits.

load(file="../Data/rData/MFHD.Rdata")

for (id in setdiff(1:15, 13)) {
  depthMFHD <- get(paste0("sub", id, "_MFHD"))$MFHDdepth
  tMFHD <- get(paste0("sub", id, "_MFHD"))$tExp
  
  cat("###", paste0("Subject", id), "{-}",'\n')
  plot(tMFHD, depthMFHD, pch=16, cex=1,col=c(rep("red", 500), rep("black", length(depthMFHD)-500)))
  legend("topright", legend=c("Baseline data", "New data"), col=c("red", "black"),
         pch=16)
  cat('\n \n')
  
}

Subject1

Subject2

Subject3

Subject4

Subject5

Subject6

Subject7

Subject8

Subject9

Subject10

Subject11

Subject12

Subject14

Subject15

4 Exponentially Weighted Moving Average (EWMA)

Cumulative SPC methods are designed to detect small and persistent changes in the data by using the observed data available at the current time point. Such methods work by accumulating the small changes in the controlled statistics that occurred prior to the present time point. Since we expected the gait data in our study to change gradually, we used the exponentially weighted moving average (EWMA) SPC method to monitor fatigue during the simulated MMH sessions.

4.1 EWMA on Subgroup Mode Depth Values

After calculating the Mode depth values for individual gait cycles we reassigned the depth values to their original subgroups in the following code chunk. Next, we used the qcc package in R and constructed the EWMA control chart by considering the \(1^{st}\) 10 subgroups as the in-control data and presented the control charts.

load(file="../Data/rData/subGs_for_depth.Rdata")
load(file="../Data/rData/mode_mag.Rdata")
aa <- 3
for (id in setdiff(1:15, 13)) {
  
  dep_Mode <- get(paste0("sub", id, "_mode_mag"))$modeDepth
  batchSize <- get(paste0("sub", id, "_aMSLT"))$batchSize
  
  ################################################ reassign depth values to original subGroups
  depMode_split <- sapply(1:length(batchSize), FUN=function(x){rep(x, batchSize[x])}) %>% do.call(what=c) %>% split(x=dep_Mode)
  
  ####################################### Create matrix with NA's at the end of smaller batches
  Mode_batch_mat <- sapply(1:length(depMode_split), FUN=function(x){c(depMode_split[[x]], rep(NA, max(batchSize) - length(depMode_split[[x]])))}) %>% t()
  
  ####################################### Remove small length batches
  rm_ind <- which(apply(Mode_batch_mat, 1, function(x) sum(!is.na(x))) < 10)
  rm_ind <- c(rm_ind, 3000)
  
  Mode_batch_mat <- Mode_batch_mat[-rm_ind,]
  
  ####################################### Plot MFHD EWMA for walking cycle batches
  Lambda <- 0.5
  end <- nrow(Mode_batch_mat)
  IC_len <- 10
  
  ewma1 <- ewma(Mode_batch_mat[1:IC_len,], lambda=Lambda, nsigmas=5, newdata=Mode_batch_mat[(IC_len+1):end,], add.stats=F, plot=F)
  
  cat("###", paste0("Subject", id), "{-}",'\n')
  plot(ewma1,
       title=paste0("subject ", id, "\nEWMA Chart on Mode Depth"), xlab="Subgroup Index")
  cat('\n \n')
}

Subject1

Subject2

Subject3

Subject4

Subject5

Subject6

Subject7

Subject8

Subject9

Subject10

Subject11

Subject12

Subject14

Subject15

4.2 EWMA on Subgroup EWMA Depth Values

In the following code chunk we have constructed the EWMA control charts on the MFHD depth values using the qcc package in R and presented the charts.

load(file="../Data/rData/subGs_for_depth.Rdata")
load(file="../Data/rData/MFHD.Rdata")

aa <- 2

for (id in setdiff(1:15, 13)) {
  
  dep_MFHD <- get(paste0("sub", id, "_MFHD"))$MFHDdepth
  batchSize <- get(paste0("sub", id, "_aMSLT"))$batchSize
  
  ################################################ reassign depth values to original subGroups
  depMFHD_split <- sapply(1:length(batchSize), FUN=function(x){rep(x, batchSize[x])}) %>% do.call(what=c) %>% split(x=dep_MFHD)
  
  ####################################### Create matrix with NA's at the end of smaller batches
  MFHD_batch_mat <- sapply(1:length(depMFHD_split), FUN=function(x){c(depMFHD_split[[x]], rep(NA, max(batchSize) - length(depMFHD_split[[x]])))}) %>% t()
  
  ####################################### Remove small length batches
  rm_ind <- which(apply(MFHD_batch_mat, 1, function(x) sum(!is.na(x))) < 10)
  rm_ind <- c(rm_ind, 3000)
  
  MFHD_batch_mat <- MFHD_batch_mat[-rm_ind,]
  
  ####################################### Plot MFHD EWMA for walking cycle batches
  Lambda <- 0.5
  end <- nrow(MFHD_batch_mat)
  IC_len <- 10
  
  ewma1 <- ewma(MFHD_batch_mat[1:IC_len,], lambda=Lambda, nsigmas=5, newdata=MFHD_batch_mat[(IC_len+1):end,], add.stats=F, plot=F)
  
  cat("###", paste0("Subject", id), "{-}",'\n')
  plot(ewma1,
       title=paste0("subject ", id, "\nEWMA Chart on MFHD Depth"), xlab="Subgroup Index")
  cat('\n \n')
}

Subject1

Subject2

Subject3

Subject4

Subject5

Subject6

Subject7

Subject8

Subject9

Subject10

Subject11

Subject12

Subject14

Subject15

---

References

Cuevas, A., Febrero, M., & Fraiman, R. (2007). Robust estimation and classification for functional data via projection-based depth notions. Computational Statistics, 22(3), 481–496.

Febrero Bande, M., & Oviedo de la Fuente, M. (2012). Statistical computing in functional data analysis: The r package fda. Usc.

Fraiman, R., & Muniz, G. (2001). Trimmed means for functional data. Test, 10(2), 419–440.

Hubert, M., Claeskens, G., De Ketelaere, B., & Vakili, K. (2012). A new depth-based approach for detecting outlying curves. Proceedings of Compstat, 329–340.

Lamooki, S. R., Kang, J., Cavuoto, L. A., Megahed, F. M., & Jones-Farmer, L. A. (n.d.). Challenges and opportunities for statistical monitoring of gait cycle acceleration observed from imu data for fatigue detection. 2020 8th Ieee Ras/Embs International Conference for Biomedical Robotics and Biomechatronics (Biorob), 593–598.

---

1. Email: saebraga@buffalo.edu | Website: LinkedIn↩︎

2. Email: jiyeonk@buffalo.edu | Phone: +1-716-645-6063 | Website: University at Buffalo Official↩︎

3. Email: loracavu@buffalo.edu | Phone: +1-716-645-4696 | Website: University at Buffalo Official↩︎

4. Email: fmegahed@miamioh.edu | Phone: +1-513-529-4185 | Website: Miami University Official↩︎

5. Email: farmerl2@miamioh.edu | Phone: +1-513-529-4823 | Website: Miami University Official↩︎